ARGUS distribution (argus)

Goal: build intuition, derive key formulas, and implement simulation/visualization for the continuous ARGUS distribution as implemented in scipy.stats.argus.

The ARGUS distribution is most famous in particle physics as a parametric model for background shapes near a hard kinematic endpoint.

import numpy as np

import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import special
from scipy.optimize import minimize_scalar
from scipy.stats import argus, chi2, norm

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

rng = np.random.default_rng(42)

print("numpy:", np.__version__)
print("scipy:", scipy.__version__)

numpy: 1.26.2
scipy: 1.15.0

1) Title & classification

• Name: argus (ARGUS distribution)

• Type: continuous

• Support (standardized): (x \in (0, 1))

  • The density is 0 at the endpoints, so people often write ([0,1]) informally.

• Parameter space: shape parameter (\chi > 0)

  • SciPy also supports loc and scale, shifting/scaling the support to ((\text{loc},, \text{loc}+\text{scale})).

2) Intuition & motivation

What it models

A typical use case is a background distribution for a quantity with a hard upper endpoint (e.g., an invariant mass that cannot exceed a known limit).

• The factor (\sqrt{1-x^2}) creates a phase-space-like suppression as (x \to 1).

• The exponential factor (\exp{-\tfrac{\chi^2}{2}(1-x^2)}) controls how strongly mass piles up near the endpoint.

Real-world use cases

• High energy physics: background modeling near kinematic limits (the classic “ARGUS function”).

• Endpoint distributions: any normalized measurement bounded above, where the density falls to 0 at the maximum.

Relations to other distributions

A very useful transformation connects ARGUS to a truncated Gamma distribution:

[ Y ;:=; \frac{\chi^2}{2}(1-X^2) \in (0, \chi^2/2). ]

Then (Y) has a density proportional to (\sqrt{y} e^{-y}), i.e. a (\text{Gamma}(3/2,,1)) conditioned on ([0,\chi^2/2]).

Limiting behavior:

• As (\chi \downarrow 0): the exponential term (\to 1), and the density approaches (f(x) \propto x\sqrt{1-x^2}).

• As (\chi \to \infty): the distribution concentrates near (x \approx 1).

3) Formal definition

PDF

For (0<x<1) and (\chi>0), the standardized ARGUS density is

[ f(x;\chi) = \frac{\chi^3}{\sqrt{2\pi},\Psi(\chi)}, x,\sqrt{1-x^2},\exp\left{-\frac{\chi^2}{2}(1-x^2)\right}, ]

where

[ \Psi(\chi) = \Phi(\chi) - \chi,\phi(\chi) - \tfrac{1}{2}, ]

and (\Phi) and (\phi) are the CDF and PDF of (\mathcal N(0,1)).

A numerically stable equivalent form is

[ \Psi(\chi) = \tfrac{1}{2},P\left(\tfrac{3}{2}, \tfrac{\chi^2}{2}\right), ]

where (P(a,x)=\gamma(a,x)/\Gamma(a)) is the regularized lower incomplete gamma function.

CDF

SciPy’s implementation uses a particularly clean survival function:

[ \bar F(x;\chi) = 1 - F(x;\chi) = \frac{\Psi\big(\chi\sqrt{1-x^2}\big)}{\Psi(\chi)}, ]

so

[ F(x;\chi) = 1 - \frac{\Psi\big(\chi\sqrt{1-x^2}\big)}{\Psi(\chi)}. ]

Location/scale

If (Z\sim\text{ARGUS}(\chi)) on ((0,1)), then (X=\text{loc}+\text{scale}\cdot Z) has support ((\text{loc},, \text{loc}+\text{scale})).

def argus_Psi(chi: np.ndarray | float) -> np.ndarray:
    """Psi(chi) used by the ARGUS distribution.

    We use the regularized incomplete gamma form for numerical stability:
        Psi(chi) = 0.5 * P(3/2, chi^2/2).
    """
    chi = np.asarray(chi, dtype=float)
    return 0.5 * special.gammainc(1.5, chi**2 / 2)


def argus_pdf(x: np.ndarray | float, chi: float) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    chi = float(chi)
    if chi <= 0:
        return np.full_like(x, np.nan, dtype=float)

    Psi = argus_Psi(chi)
    norm_const = chi**3 / (np.sqrt(2 * np.pi) * Psi)

    y = 1.0 - x * x
    base = norm_const * x * np.sqrt(np.clip(y, 0.0, None)) * np.exp(-0.5 * chi**2 * y)
    return np.where((x > 0) & (x < 1), base, 0.0)


def argus_cdf(x: np.ndarray | float, chi: float) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    chi = float(chi)
    if chi <= 0:
        return np.full_like(x, np.nan, dtype=float)

    Psi_chi = argus_Psi(chi)
    t = chi * np.sqrt(np.clip(1.0 - x * x, 0.0, None))
    sf = argus_Psi(t) / Psi_chi

    cdf = 1.0 - sf
    cdf = np.where(x <= 0, 0.0, cdf)
    cdf = np.where(x >= 1, 1.0, cdf)
    return cdf


# Quick consistency check vs SciPy
xgrid = np.linspace(0, 1, 400)
chi_test = 2.5

max_pdf_err = np.max(np.abs(argus_pdf(xgrid, chi_test) - argus.pdf(xgrid, chi_test)))
max_cdf_err = np.max(np.abs(argus_cdf(xgrid, chi_test) - argus.cdf(xgrid, chi_test)))

print("max |pdf - scipy|:", max_pdf_err)
print("max |cdf - scipy|:", max_cdf_err)

max |pdf - scipy|: 1.7763568394002505e-15
max |cdf - scipy|: 3.552713678800501e-15

4) Moments & properties

Mean and variance

SciPy uses (and we can derive) the following closed forms.

Let (I_1) be the modified Bessel function of the first kind. Define (z=\chi^2/4).

[ \mathbb E[X] = \sqrt{\frac{\pi}{8}},\frac{\chi,e^{-z} I_1(z)}{\Psi(\chi)}. ]

For the second moment: [ \mathbb E[X^2] = 1 - \frac{3}{\chi^2} + \frac{\chi,\phi(\chi)}{\Psi(\chi)}. ]

Then (\operatorname{Var}(X)=\mathbb E[X^2]-\mathbb E[X]^2).

Skewness and kurtosis

Skewness and (excess) kurtosis are defined via central moments: [ \gamma_1 = \frac{\mu_3}{\sigma^3},\qquad \gamma_2 = \frac{\mu_4}{\sigma^4} - 3. ]

For ARGUS, these don’t simplify nicely to a short expression; they’re typically computed numerically.

MGF / characteristic function

Because (X\in(0,1)) is bounded, both exist for all real arguments: [ M_X(t)=\mathbb E[e^{tX}],\qquad \varphi_X(\omega)=\mathbb E[e^{i\omega X}]. ]

There is no widely used simple closed form; numerical quadrature (or a moment series) is standard.

Entropy

The differential entropy is [ H(X) = -\int_0^1 f(x;\chi),\log f(x;\chi),dx. ]

SciPy provides argus.entropy(chi).

def argus_mean(chi: np.ndarray | float) -> np.ndarray:
    chi = np.asarray(chi, dtype=float)
    Psi = argus_Psi(chi)
    z = chi**2 / 4
    return np.sqrt(np.pi / 8) * chi * special.ive(1, z) / Psi


def argus_E_x2(chi: np.ndarray | float) -> np.ndarray:
    """E[X^2] with a small-chi branch for numerical stability (mirrors SciPy)."""
    chi = np.asarray(chi, dtype=float)
    out = np.empty_like(chi)

    mask = chi > 0.1
    if np.any(mask):
        c = chi[mask]
        out[mask] = 1 - 3 / c**2 + c * norm.pdf(c) / argus_Psi(c)

    if np.any(~mask):
        c = chi[~mask]
        # series approximation for small chi (from SciPy's implementation)
        coef = [-358 / 65690625, 0, -94 / 1010625, 0, 2 / 2625, 0, 6 / 175, 0, 0.4]
        out[~mask] = np.polyval(coef, c)

    return out


def argus_var(chi: np.ndarray | float) -> np.ndarray:
    m = argus_mean(chi)
    return argus_E_x2(chi) - m**2


chis = np.array([0.2, 0.5, 1.0, 2.5, 6.0])

m_formula = argus_mean(chis)
v_formula = argus_var(chis)

m_scipy, v_scipy, s_scipy, k_scipy = argus.stats(chis, moments="mvsk")

print("chi    mean(formula)   mean(scipy)   var(formula)    var(scipy)   skew     kurt(excess)")
for chi, mf, ms, vf, vs, ss, ks in zip(chis, m_formula, m_scipy, v_formula, v_scipy, s_scipy, k_scipy):
    print(f"{chi:4.1f}  {mf:12.6f}  {ms:11.6f}  {vf:12.6f}  {vs:11.6f}  {ss:7.4f}  {ks:11.4f}")


# MGF/CF: demonstrate numerical quadrature vs Monte Carlo
chi0 = 2.5
x_mc = argus.rvs(chi0, size=150_000, random_state=rng)

ts = np.array([-3.0, 0.0, 3.0])
ws = np.array([0.0, 10.0, 20.0])

print("\nMGF M(t)=E[e^{tX}] at a few t:")
for t in ts:
    mc = np.mean(np.exp(t * x_mc))
    quad = argus.expect(lambda x, t=t: np.exp(t * x), args=(chi0,))
    print(f"t={t:>5.1f}  MC={mc:.6f}  quad={quad:.6f}")

print("\nCharacteristic function φ(ω)=E[e^{iωX}] at a few ω:")
for w in ws:
    mc = np.mean(np.exp(1j * w * x_mc))
    quad = argus.expect(lambda x, w=w: np.exp(1j * w * x), args=(chi0,))
    print(f"w={w:>5.1f}  MC={mc:.6f}  quad={quad:.6f}")

print("\nEntropy at chi=2.5:", argus.entropy(chi0))

chi    mean(formula)   mean(scipy)   var(formula)    var(scipy)   skew     kurt(excess)
 0.2      0.590227     0.590227      0.053005     0.053005  -0.2966      -0.7982
 0.5      0.596428     0.596428      0.052891     0.052891  -0.3230      -0.7800
 1.0      0.618703     0.618703      0.052157     0.052157  -0.4204      -0.6939
 2.5      0.762779     0.762779      0.035555     0.035555  -1.1882       1.0606
 6.0      0.956717     0.956717      0.001359     0.001359  -1.8768       5.8944

MGF M(t)=E[e^{tX}] at a few t:
t= -3.0  MC=0.123343  quad=0.123745
t=  0.0  MC=1.000000  quad=1.000000
t=  3.0  MC=11.224313  quad=11.216045

Characteristic function φ(ω)=E[e^{iωX}] at a few ω:
w=  0.0  MC=1.000000+0.000000j  quad=1.000000
w= 10.0  MC=-0.252645+0.292744j  quad=-0.253406
w= 20.0  MC=0.133822-0.107030j  quad=0.133083

Entropy at chi=2.5: -0.48713485149744684

/home/tempa/miniconda3/lib/python3.12/site-packages/scipy/integrate/_quadpack_py.py:606: ComplexWarning:

Casting complex values to real discards the imaginary part

5) Parameter interpretation (shape changes)

The single shape parameter (\chi) controls how quickly the density falls away from the endpoint.

• Small (\chi): the exponential term is weak, and the shape is dominated by (x\sqrt{1-x^2}).

• Large (\chi): the exponential term strongly favors (x) near 1, creating a sharp peak close to the endpoint.

A convenient closed-form for the mode (maximum of the pdf) comes from differentiating (\log f(x)):

[ \text{mode}(X)^2 = \frac{\chi^2 - 2 + \sqrt{\chi^4 + 4}}{2\chi^2}. ]

def argus_mode(chi: float) -> float:
    chi = float(chi)
    if chi <= 0:
        return np.nan
    u = (chi**2 - 2 + np.sqrt(chi**4 + 4)) / (2 * chi**2)
    return float(np.sqrt(u))


for chi in [0.2, 1.0, 2.5, 6.0]:
    print(f"chi={chi:>4}: mode≈{argus_mode(chi):.4f}")

chi= 0.2: mode≈0.7106
chi= 1.0: mode≈0.7862
chi= 2.5: mode≈0.9300
chi= 6.0: mode≈0.9864

6) Derivations

6.1 Derivation of the mean (why Bessel functions appear)

Start from ( \mathbb E[X] = \int_0^1 x,f(x;\chi),dx )

so the integral of interest is

[ I = \int_0^1 x^2\sqrt{1-x^2},\exp\left{-\frac{\chi^2}{2}(1-x^2)\right},dx. ]

Use the substitution (x=\cos\theta) with (\theta\in[0,\pi/2]). Then (1-x^2=\sin^2\theta), (dx=-\sin\theta,d\theta), and

[ I = \int_0^{\pi/2} \cos^2\theta,\sin^2\theta,\exp\left{-\frac{\chi^2}{2}\sin^2\theta\right}d\theta. ]

Rewrite (\sin^2\theta = \tfrac{1-\cos 2\theta}{2}) and set (\varphi=2\theta), which yields an integral of the form (\int_0^{\pi} e^{z\cos\varphi}\cos(n\varphi)d\varphi), whose value is (\pi I_n(z)).

After simplification you arrive at

[ I = \frac{\pi}{2\chi^2} e^{-\chi^2/4} I_1(\chi^2/4). ]

Multiplying by the normalization constant (\chi^3/(\sqrt{2\pi}\Psi(\chi))) gives

[ \mathbb E[X] = \sqrt{\frac{\pi}{8}},\frac{\chi,e^{-\chi^2/4}I_1(\chi^2/4)}{\Psi(\chi)}. ]

6.2 Derivation of (\mathbb E[X^2]) via the truncated-Gamma representation

Let (Y=\tfrac{\chi^2}{2}(1-X^2)). Then (Y) is Gamma((3/2,1)) conditioned on ([0,\chi^2/2]).

Since (X^2=1-2Y/\chi^2):

[ \mathbb E[X^2] = 1 - \frac{2}{\chi^2},\mathbb E[Y\mid Y\le \chi^2/2]. ]

For a Gamma((\alpha,1)) truncated at (b), (\mathbb E[Y\mid Y\le b] = \gamma(\alpha+1,b)/\gamma(\alpha,b)).

Using the incomplete-gamma recurrence simplifies the ratio and leads to

[ \mathbb E[X^2] = 1 - \frac{3}{\chi^2} + \frac{\chi,\phi(\chi)}{\Psi(\chi)}. ]

6.3 Likelihood for i.i.d. data

For data (x_1,\dots,x_n\in(0,1)), the log-likelihood (standardized) is

[ \ell(\chi) = n\left(3\log\chi - \log\Psi(\chi) - \tfrac{1}{2}\log(2\pi)\right)

• \sum_{i=1}^n\left(\log x_i + \tfrac{1}{2}\log(1-x_i^2) - \tfrac{\chi^2}{2}(1-x_i^2)\right). ]

A handy derivative identity is (\Psi’(\chi)=\chi^2\phi(\chi)), so you can optimize (\ell(\chi)) with gradient-based methods.

def argus_loglik(chi: float, x: np.ndarray) -> float:
    """Log-likelihood for standardized ARGUS(chi) given x in (0,1)."""
    chi = float(chi)
    if chi <= 0:
        return -np.inf

    x = np.asarray(x, dtype=float)
    if np.any((x <= 0) | (x >= 1)):
        return -np.inf

    Psi = argus_Psi(chi)

    y = 1.0 - x * x
    ll = (
        x.size * (3 * np.log(chi) - 0.5 * np.log(2 * np.pi) - np.log(Psi))
        + np.sum(np.log(x) + 0.5 * np.log1p(-x * x) - 0.5 * chi**2 * y)
    )
    return float(ll)


def argus_loglik_grad(chi: float, x: np.ndarray) -> float:
    """Gradient of the log-likelihood.

    Uses Psi'(chi) = chi^2 * phi(chi).
    """
    chi = float(chi)
    if chi <= 0:
        return np.nan

    x = np.asarray(x, dtype=float)
    if np.any((x <= 0) | (x >= 1)):
        return np.nan

    n = x.size
    Psi = argus_Psi(chi)
    Psi_prime = chi**2 * norm.pdf(chi)

    y = 1.0 - x * x
    return float(n * (3 / chi - Psi_prime / Psi) - chi * np.sum(y))


# Simulate data and compute the MLE
chi_true = 2.5
x = argus.rvs(chi_true, size=800, random_state=rng)

obj = lambda c: -argus_loglik(c, x)
res = minimize_scalar(obj, bounds=(1e-6, 30.0), method="bounded")
chi_hat = float(res.x)

print("chi_true:", chi_true)
print("chi_hat (MLE):", chi_hat)
print("grad at chi_hat:", argus_loglik_grad(chi_hat, x))

# Likelihood ratio test: H0: chi = chi0 vs H1: chi free
chi0 = 1.0
lr_stat = 2 * (argus_loglik(chi_hat, x) - argus_loglik(chi0, x))
p_value = chi2.sf(lr_stat, df=1)

print("\nLRT statistic:", lr_stat)
print("approx p-value (chi-square df=1):", p_value)

chi_true: 2.5
chi_hat (MLE): 2.5179774750503756
grad at chi_hat: -2.3556725523121713e-05

LRT statistic: 386.62155391177384
approx p-value (chi-square df=1): 4.5017064542421926e-86

7) Sampling & simulation (NumPy-only)

Below is a NumPy-only sampler for standardized ARGUS((\chi)). It mirrors SciPy’s piecewise strategy:

• Small (\chi): rejection sampling using the (\chi\to 0) base density (g(x)=3x\sqrt{1-x^2}).

• Moderate (\chi): rejection sampling with a proposal density (g(x)\propto x\exp{-\tfrac{\chi^2}{2}(1-x^2)}).

• Large (\chi): use the truncated Gamma representation of (Y=\tfrac{\chi^2}{2}(1-X^2)).

This is not the only way to sample ARGUS, but it’s easy to implement and performs well across a wide parameter range.

def argus_rvs_numpy(chi: float, size=1, rng: np.random.Generator | None = None) -> np.ndarray:
    """Sample from standardized ARGUS(chi) using only NumPy."""
    chi = float(chi)
    if chi <= 0:
        raise ValueError("chi must be > 0")

    if rng is None:
        rng = np.random.default_rng()

    size1d = tuple(np.atleast_1d(size))
    n = int(np.prod(size1d))

    out = np.empty(n, dtype=float)
    simulated = 0

    chi2 = chi * chi

    if chi <= 0.5:
        # Case 1: propose from g(x) = 3*x*sqrt(1-x^2), accept with exp(-chi^2(1-x^2)/2)
        d = -chi2 / 2
        while simulated < n:
            k = n - simulated
            u = rng.uniform(size=k)
            v = rng.uniform(size=k)
            z = v ** (2 / 3)  # z = 1 - x^2 under the proposal
            accept = np.log(u) <= d * z
            num_accept = int(np.sum(accept))
            if num_accept:
                out[simulated : simulated + num_accept] = np.sqrt(1 - z[accept])
                simulated += num_accept

    elif chi <= 1.8:
        # Case 2: propose from g(x) ∝ x*exp(-chi^2(1-x^2)/2) and accept with sqrt(1-x^2)
        echi = np.exp(-chi2 / 2)
        while simulated < n:
            k = n - simulated
            u = rng.uniform(size=k)
            v = rng.uniform(size=k)

            # z <= 0, and x = sqrt(1 + z)
            z = 2 * np.log(echi * (1 - v) + v) / chi2
            accept = (u * u + z) <= 0
            num_accept = int(np.sum(accept))
            if num_accept:
                out[simulated : simulated + num_accept] = np.sqrt(1 + z[accept])
                simulated += num_accept

    else:
        # Case 3: conditional Gamma(3/2, 1) for Y in [0, chi^2/2]
        y = np.empty(n, dtype=float)
        while simulated < n:
            k = n - simulated
            g = rng.standard_gamma(shape=1.5, size=k)  # Gamma(k=3/2, theta=1)
            accept = g <= chi2 / 2
            num_accept = int(np.sum(accept))
            if num_accept:
                y[simulated : simulated + num_accept] = g[accept]
                simulated += num_accept
        out = np.sqrt(1 - 2 * y / chi2)

    return out.reshape(size1d)


# Sanity check: NumPy sampler vs SciPy moments
chi_check = 2.5
s = argus_rvs_numpy(chi_check, size=200_000, rng=rng)

print("sample mean:", s.mean())
print("theory mean:", argus.mean(chi_check))

print("sample var:", s.var())
print("theory var:", argus.var(chi_check))

sample mean: 0.7621546176726476
theory mean: 0.7627785650581257
sample var: 0.035698472048319115
theory var: 0.035554890425441354

8) Visualization (PDF, CDF, Monte Carlo)

We’ll visualize how (\chi) changes the shape and verify the sampler by overlaying a Monte Carlo histogram on the theoretical PDF.

x = np.linspace(0, 1, 800)
chis_plot = [0.2, 0.7, 2.5, 6.0]

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=("PDF", "CDF"),
)

for chi in chis_plot:
    fig.add_trace(go.Scatter(x=x, y=argus_pdf(x, chi), name=f"chi={chi}", mode="lines"), row=1, col=1)
    fig.add_trace(go.Scatter(x=x, y=argus_cdf(x, chi), name=f"chi={chi}", mode="lines", showlegend=False), row=1, col=2)

fig.update_xaxes(title_text="x (standardized)", row=1, col=1)
fig.update_xaxes(title_text="x (standardized)", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="CDF", row=1, col=2)
fig.update_layout(title="ARGUS distribution for different chi", width=950)
fig.show()


# Monte Carlo overlay
chi_mc = 2.5
n_mc = 40_000
samples = argus_rvs_numpy(chi_mc, size=n_mc, rng=rng)

hist = np.histogram(samples, bins=70, range=(0, 1), density=True)
bins = hist[1]
centers = 0.5 * (bins[:-1] + bins[1:])

fig2 = go.Figure()
fig2.add_trace(go.Bar(x=centers, y=hist[0], name="MC histogram", opacity=0.4))
fig2.add_trace(go.Scatter(x=x, y=argus_pdf(x, chi_mc), name="theoretical pdf", mode="lines"))
fig2.update_layout(
    title=f"Monte Carlo check (chi={chi_mc})",
    xaxis_title="x",
    yaxis_title="density",
    bargap=0.02,
)
fig2.show()

9) SciPy integration (scipy.stats.argus)

SciPy exposes ARGUS as a standard rv_continuous distribution:

• argus.pdf(x, chi) / argus.logpdf(x, chi)

• argus.cdf(x, chi) / argus.sf(x, chi)

• argus.rvs(chi, size=..., random_state=...)

• argus.fit(data, ...)

Because the support is bounded, it’s common to pre-normalize data to ((0,1)) and fit only (\chi) by fixing loc=0, scale=1.

chi = 2.5
x = np.linspace(0, 1, 6)
print("x:", x)
print("pdf:", argus.pdf(x, chi))
print("cdf:", argus.cdf(x, chi))

# Random variates
r = argus.rvs(chi, size=5, random_state=rng)
print("\nrvs:", r)

# Fit chi (fix loc/scale)
chi_true = 2.5
x_data = argus.rvs(chi_true, size=2000, random_state=rng)

chi_hat, loc_hat, scale_hat = argus.fit(x_data, floc=0.0, fscale=1.0)
print("\ntrue chi:", chi_true)
print("fit chi :", chi_hat)
print("(loc, scale fixed to)", loc_hat, scale_hat)

x: [0.  0.2 0.4 0.6 0.8 1. ]
pdf: [0.         0.13515406 0.36789472 0.89991027 2.15877251 0.        ]
cdf: [0.         0.01283353 0.06035862 0.17934912 0.46903877 1.        ]

rvs: [0.44379125 0.96750725 0.75501664 0.64295312 0.792296  ]

true chi: 2.5
fit chi : 2.528417968750004
(loc, scale fixed to) 0.0 1.0

10) Statistical use cases

10.1 Hypothesis testing

A common workflow is to compare a fixed-shape ARGUS background (e.g. (\chi=\chi_0)) against a fitted (\chi) using a likelihood ratio test.

Caveat: the usual (\chi^2) calibration is asymptotic and can be inaccurate for small samples.

10.2 Bayesian modeling

Treat (\chi) as an unknown parameter with a prior (e.g. log-normal). Because it’s 1D, a grid posterior is often sufficient.

10.3 Generative modeling

ARGUS is frequently used as a background component inside a mixture model:

[ p(x)=\pi,p_{\text{signal}}(x) + (1-\pi),p_{\text{bkg}}(x),\qquad p_{\text{bkg}}(x)=\text{ARGUS}(\chi). ]

Below is a lightweight demo of all three.

# --- 10.1 Likelihood ratio test demo ---
chi_true = 2.5
x = argus.rvs(chi_true, size=600, random_state=rng)

res = minimize_scalar(lambda c: -argus_loglik(c, x), bounds=(1e-6, 30.0), method="bounded")
chi_hat = float(res.x)

chi0 = 1.0
lr_stat = 2 * (argus_loglik(chi_hat, x) - argus_loglik(chi0, x))
print("chi_hat:", chi_hat)
print("LRT stat:", lr_stat)
print("approx p-value:", chi2.sf(lr_stat, df=1))


# --- 10.2 Bayesian grid posterior for chi ---
# Prior: log chi ~ Normal(mu=0, sigma=1)  (=> chi is log-normal)
chis = np.linspace(0.05, 10.0, 400)
log_prior = norm.logpdf(np.log(chis), loc=0.0, scale=1.0) - np.log(chis)  # Jacobian for chi -> log chi
log_like = np.array([argus_loglik(c, x) for c in chis])

log_post_unnorm = log_like + log_prior
log_post_unnorm -= np.max(log_post_unnorm)
post = np.exp(log_post_unnorm)
post /= np.trapz(post, chis)

post_cdf = np.cumsum(post)
post_cdf /= post_cdf[-1]

def quantile(q):
    return float(np.interp(q, post_cdf, chis))

ci_low, ci_high = quantile(0.05), quantile(0.95)
post_mean = float(np.trapz(chis * post, chis))

print("\nPosterior mean:", post_mean)
print("90% credible interval:", (ci_low, ci_high))

fig = go.Figure(go.Scatter(x=chis, y=post, mode="lines", name="posterior"))
fig.add_vline(x=chi_true, line_dash="dash", line_color="black", annotation_text="true")
fig.add_vrect(x0=ci_low, x1=ci_high, fillcolor="lightblue", opacity=0.3, line_width=0)
fig.update_layout(title="Posterior over chi (grid)", xaxis_title="chi", yaxis_title="density")
fig.show()


# --- 10.3 Simple generative mixture (signal + ARGUS background) ---
# Background: ARGUS
chi_bkg = 6.0
n = 50_000

# Signal: a narrow truncated normal near the endpoint
mu_sig, sigma_sig = 0.97, 0.02
pi_sig = 0.08

n_sig = int(round(pi_sig * n))
n_bkg = n - n_sig

x_bkg = argus_rvs_numpy(chi_bkg, size=n_bkg, rng=rng)

# Truncated normal via rejection (fine for a demo)
x_sig = []
while len(x_sig) < n_sig:
    z = rng.normal(loc=mu_sig, scale=sigma_sig, size=n_sig)
    z = z[(z > 0) & (z < 1)]
    x_sig.extend(z.tolist())
x_sig = np.array(x_sig[:n_sig])

x_mix = np.concatenate([x_bkg, x_sig])

fig = px.histogram(
    x_mix,
    nbins=120,
    histnorm="probability density",
    title="Mixture example: ARGUS background + truncated-normal signal",
    labels={"value": "x"},
)
xx = np.linspace(0, 1, 800)
fig.add_trace(
    go.Scatter(
        x=xx,
        y=(1 - pi_sig) * argus_pdf(xx, chi_bkg),
        name="(1-π) * ARGUS pdf",
        mode="lines",
    )
)
fig.update_layout(showlegend=True)
fig.show()

chi_hat: 2.526546565812786
LRT stat: 294.2642747018417
approx p-value: 5.853726097616042e-66

Posterior mean: 2.5215534518729363
90% credible interval: (2.3998870486119035, 2.6156776669782578)

11) Pitfalls

• Parameter constraints: (\chi\le 0) is invalid.

• Boundary values: the theoretical support is ((0,1)). Real data may contain exact 0/1 due to rounding; for likelihood work you may need to clip slightly, e.g. x = np.clip(x, 1e-12, 1-1e-12).

• Numerical stability near 1: use log1p(-x*x) rather than log(1-x**2).

• Sampling efficiency: the simple truncated-Gamma method is great for large (\chi) but can reject a lot when (\chi) is small; the piecewise sampler avoids that.

• Fitting (fit) surprises: by default, SciPy also fits loc and scale. If your data are already standardized, fix them with floc=0, fscale=1.

12) Summary

• argus is a continuous distribution on ((0,1)) with shape parameter (\chi>0).

• Its PDF combines a phase-space term (x\sqrt{1-x^2}) and an exponential tilt controlled by (\chi).

• The CDF is especially clean in terms of (\Psi(\chi)): (F(x)=1-\Psi(\chi\sqrt{1-x^2})/\Psi(\chi)).

• Mean and second moment have closed forms (Bessel/incomplete-gamma), while skewness/kurtosis are usually computed numerically.

• Sampling can be implemented with rejection methods and a truncated-Gamma transformation; SciPy provides a fast and robust implementation in scipy.stats.argus.